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Abstract 

We present precision Monte Carlo calculations solving the QCD evolution equa- 
tions up to the next-to-leading-order (NLO) level. They employ forward Markovian 
Monte Carlo (FMC) algorithms, which provide the rigorous solutions of the QCD 
evolution equations. Appropriate Monte Carlo algorithms are described in detail. 
They are implemented in the form of the Monte Carlo program EvolFMC, which 
features the NLO kernels for the QCD evolution. The presented numerical results 
agree with those from independent, non-MC, programs (QCDNuml6, APCheb33) at 
the level of 0.1%. In this way we have demonstrated the feasibility of the precision 
MC calculations for the QCD evolution and provided very useful numerical tests 
(benchmarks) for other, non-Markovian, MC algorithms developed recently. 
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1 Introduction 



It is commonly known that the so-called evolution equations of the quark and gluon 
distributions in the hadron, derived in QED and QCD using the renormalization group or 
diagrammatic techniques [1], can be interpreted probabilistically as a Markovian process, 
see e.g. Ref. [2]. Such a process can be modeled using Monte Carlo (MC) methods. The 
corresponding MC algorithm, called in the following the Markovian MC, provides, in 
principle, an exact solution of the evolution equations for parton distribution functions 
(PDFs). In practice, the main limitation of a such solution is the size of a generated 
MC sample, i.e. corresponding statistical errors of numerical results. This is probably the 
main reason why this possibility has not been exploited until recently. Instead, alternative 
numerical methods and programs solving the QCD evolution equations much faster than 
the Markovian MC have been used. Typical examples of such non-MC programs are 
QCDNumie [3] and APCheb33 [4], see also Ref. [5]. 

Feasibility of solving efficiently the DGLAP equations [1] at the leading-order (LO) 
approximation with the Markovian MC was demonstrated for the first time by two of 
us (S.J. and M.S.) in Ref. [6]. There, the basic formahsm was briefly sketched and 
the first numerical results were presented. Good agreement between the constructed 
Markovian MC program and QCDNuinl6 for gluon and quark-singlet distribution functions 
was achieved. However, some small residual differences, at the level of 0.2%, between 
the two programs were found. Their origin was not understood at that time. Here we 
repeat the above comparisons, explain the source of these discrepancies and show the 
corrected results which agree at the level of 0.1%. The main conclusion of Ref. [6] was 
that the currently available computer CPU power allows to solve efficiently and precisely 
(at the per-mill level) the QCD evolution equations with the use of the Markovian MC 
algorithm. Of course, this method will always be slower in CPU time than non-MC 
techniques. However, it has several advantages, such as: no biases and/or numerical 
instabilities related to finite grids of points, use of quadratures, decomposition into finite 
series of polynomials, accumulation of rounding errors, etc. It is also more fiexible in 
treatment of PDFs (e.g. no need to split them into singlet and non-singlet components) 
and easier to extend into higher orders, new contributions, etc. 

The above Markovian algorithm can be a basis for the final-state radiation (FSR) 
parton shower MC program that not only solves numerically the evolution equations 
but also generates events in terms of parton flavours and four-momenta. Moreover, this 
algorithm can be a starting point and a testing tool for various kinds of constrained MC 
algorithms [7-10] being developed for the initial-state radiation (ISR). 

This paper is devoted to the Markovian MC solution of the DGLAP evolution equa- 
tions up to the next-to-leading order (NLO) in the perturbative QCD. It is organized as 
follows: In Section 2 we present a general structure of the DGLAP equations and dis- 
cuss their basic features up to the next-to-next-to- leading order (NNLO). In Section 3 
we describe in detail the Markovian MC algorithm for parton density distributions. We 
start from a classic iterative solution of the DGLAP equations and show how it can be 
expressed in terms of Markovian transition probabilities. Then we provide a method for 
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generation of a single Markovian step according to these probabilities. Subsection 3.4 is 
devoted to construction of a weighted Markovian MC algorithm, where some importance 
sampling is used to generate evolution variables. In the last subsection we show how 
the above algorithm can be modified in order to account for the running QCD coupling 
constant. In Section 4 we present the Markovian MC algorithm for parton- momentum 
distributions. It has certain advantages over the previous algorithm due to momentum 
sum rules that can be applied to evolution kernels. Both the above algorithms have been 
implemented in the MC program called EvolFMC [11]. Numerical results from EvolFMC at 
the LO and the NLO are presented in Section 5. They are compared with the results of 
non-MC programs QCDNuml6 and APCheb33. Section 6 summarizes the paper and gives 
some outlook for the future. In Appendix A we collect formulae for the QCD kernels 
(splitting functions) up to the NLO as well as explicit expressions for the NLO Sudakov 
form-factor. Appendices B and C contain formulae for simplified evolution kernels that 
are used for importance sampling in the weighted Markovian algorithm. Finally, in Ap- 
pendix D we discuss a generic discrete Markovian process. It can be seen as an illustration 
of basic features of the DGLAP-like evolution equations and their solution in terms of the 
Markovian MC algorithm. 

2 QCD evolution equations 

2.1 General structure of DGLAP equations 

The DGLAP evolution equations for quark, antiquark and gluon distributions, 

{?!,..., gi,...,g„^, G}(ii,x), (1) 
take the following general form [1, 12] 



d 

Qi 



9 In fj? 



9 In /x^ ' 
d 



G 



Qi = Y1 (^5,., ® Qj + P-m, ® Qj) + P-,,G 
j 

Yl {^Gg, (8) Qj + Pa-q, ® Qj) + Pgg^G (2) 



9 In fjL 



G 



where the summation is performed over quark flavours, j = 1, . . . .nf. The parton dis- 
tributions are functions of the Bjorken variable x and the factorization scale ^, identified 
with a hard scale in a given process (e.g. n = ^/Q^ in deep inelastic scattering). The func- 
tions P = P{ii,x) are splitting functions to be discussed below. The integral convolution 
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denoted by ® involves only longitudinal momentum fractions 

1 1 

{P ig) q){fi,x) = j dy j dzS{x - zy) P{as,z) q{fi,y) 



1 1 

/dz / X \ f dz / X \ 

— P{as,z)qi^^,-j = j — P [as, -j q{li, z) . (3) 

X X 

The splitting functions P{as, z) depend on /i through the strong coupling constant = 

P(a„.) = gpW(.) + {^yP''\z) + {^yP''\z) +.... (4) 

The superscripts (0), (1), (2) refer respectively to the leading (LO), next-to-leading (NLO) 
and next-to-next-to-leading order (NNLO) approximations in which the splitting functions 
are computed^. 

From charge conjugation and SU{nf) symmetry the splitting functions P have the 
following general structure which is independent of the approximation in which they have 
been computed 



p 


= P-- 


"JJ-* qq 


+ P^ 

~ qq 


P _ 


— P_ 




+ P^- 
~ qq 


Pq^G 


= Pq^G 


= PpG 




Pcq, 


= PGq, 


= Pgf ■ 





(5) 

Substituting these relations to Q, we find 
d 



Qi = Pqq ®qi + Pqq ® Qi + Pqq ® + ® + PpG ® G 

3 3 

ain/i^^^ = Pqq®^^ + Pqq®^^ + Pq-q®Y.'i^+Pqq®Y.^3+PFG®G 



9 In /i^ 

d 



3 

d 



d\na^^ = PGF(^Y.^q,+q,) + Pgg(^G (6) 

3 

This is the basic form of the DGLAP evolution equations. 

Within a given approximation some splitting functions may vanish or be equal. In 
particular, 



"'^We adopt the convention of Curei, Furmanski and Petronzio [13, 14] in which the expansion parameter 
equals as/(27r). The NNLO analysis of Moch and Vogt [15, 16] uses as/ {An). 
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in the LO [1] 

pV{0) ^ p5(0) ^ pS{0) ^ p 

^qq qq qq ^ ' \' ) 

in the NLO [13, 14] 



but in the NNLO [15, 16] 

qq ^ qq 



PqT ^ PS^- (9) 



Eqs. © can be rewritten in an alternative form which involves quark singlet and 
non-singlet distributions. We will present this form below. 

2.1.1 Singlet case 

The quark singlet distribution is defined as 

Uf 

= ■ (10) 

Performing summation over quark flavours in the first two equations ©, we find 

S = {Pj^ + Pj^ + nf{P^^ + Pf,)} ® S + {2nfP,c) ® G (11) 
Introducing the notation 

PFF = PX + nfPl (12) 

py/ = pyf + py.^^ (13) 

the following closed set of equations is obtained for the quark singlet and gluon distribu- 
tions 

-^-^ ^ = Pff®^ + i^UfPFo) ® G (14) 



(9 In /i 

d 



G = Pgf®^ + Pgg®G. (15) 



(9 In /i^ 

The splitting functions in these equations obey the general relations 

1 1 
j dz{zPFF{l^,z) + zPgf{iJ',z)} = J dz{2nfzPFGifJ',z) + zPGGilJ',z)} = 0. (16) 



4 



They immediately leads to the momentum sum rule 



dx {x'Z{fi,x) + xG{fi,x)} = const. (17) 



which is conserved during the evolution. In the parton model interpretation the constant 
is set to one by normalizing the initial conditions for Eqs. (jl4ll5|) : T,{fio,x) and G{fio,x). 

2.1.2 Non-singlet case 

Introducing the quark non-singlet distribution 

y{^^^x) = -qj{fi,x)) , (18) 

i=i 

the following evolution equation is obtained from Eqs. © 

where the new splitting function reads 

P^s = P-+ njP'l (20) 

pyj = pyf-py^s_ (21) 

Similarly, for the non-singlet quark distributions 

Qii^^x) = -Qiif^,^) - —V{fi,x) (22) 

^/ 

iHfJ^^ ^) = 9i(/^, x) + qi{^, x) - — S(/i, x) , (23) 

we find from Eqs. (fT^ and (fT^ the following equations 

d 



d 



q- = py_®q- (24) 
-qf = Pl® qf . (25) 



9 In yU- 

Notice that there is no gluon distribution in the derived equations 
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2.2 Summary of the forms 

With the sphtting functions usually presented in the literature 

{Pjj -Pf ' Pfg, Pgf, Pgg} , (26) 



the evolution equations for the parton distributions {q^ , q^, V, S, G} are given by Eqs. 
(123), (Cni), (El) and (US)), respectively. 

According to relations (I7j)-(jni), the perturbative expansions for the kernels and 
take the following form 



P+ [as, z 



a. 
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P'-{as,z)= P'}'\z) + ... (27) 

The remaining kernels {Pfg, Pgf, Pgg} have the nonzero splitting functions in each 
approximation. 

Alternatively, the parton distributions {qi, q^, G} could be evolved with the help of 
Eqs. (jni) with the kernels 

{Pqq^y Pqi^-> Pfg, Pgf, Pgg} , (28) 
where P^^^ and P^^ are computed by inverting relations ()13p and ()2ip: 

pv,s ^ i ^pv,s _ pv,s^^ (30) 
2.3 Behaviour at z ^ 1 

Let us consider the splitting functions (j26j) . All the kernels, except the Pf, are divergent 
for z = 1. 

The quark-quark and gluon-gluon sphtting functions {P±, P^, Pgg} have the follow- 
ing form 

P(a„ z) = -^^^ + Bias) 5(1 - z) + P(a„ z) , (31) 
(1-2;)+ 

where the "+" prescription regularizes the 1/(1 — 2;) singularity 

1 

[f{z)U = f{z)-S{l-z) fdz'fiz')., (32) 



and the functions A{as), B{as) and P{as,z) are computed in powers of see Eq. (jH). 
In particular 

P(a„z) = J^^s'-'D^'^z). (33) 

fc=0 

In the LO approximation (k = 0) D^^\z = 1) is finite [1] while in the NLO {k = 1) and 
NNLO {k = 2) approximations, the coefficients are logarithmically divergent [13-16]: 



D^''\z) = Dk In(l-z) + 0{1). 



(34) 



Similarly, the quark-gluon and gluon-quark splitting functions {Pfg, Pgf} contain 
logarithmically divergent terms for z = 1 [14, 16]: 



fe=0 



2k 



i=l 



(35) 



P{a„z) 



(36) 



Thus in the limit z — 1, we have for Ppc and Pgf'- 

' 0{as) in LO (k = 0) 

0{a'i ln\l - z)) inNLO(k=l) 

^ 0{alln\l - z)) in NNLO (k = 2). 

2.4 Behaviour at 2; ^ 

As in the previous section, let us consider the splitting functions 

The splitting functions {P^, P^} are logarithmically divergent at z = starting from 
the NLO approximation: [13-16]: 



P(a.„ z) 



E"^M5^A^'' ln^^ + 0(l)}, 

k=0 ^ i=l 



(37) 



Thus for 2; — > 0, we find for P^ and P^: 



0{a^Jn^z) 
^ 0{a^Jn^z) 



P(a., z) 



in LO (k = 0) 
in NLO (k = 1) 
in NNLO (k = 2). 



(3^ 



The remaining splitting functions {P^,Pfg, Pgf, Pgg} have the following behaviour 
for 2 ^ [14,16]: 



P{as,z) = Ei{a, 



^^hl + E,{a,)- + 0{\n''z) 

z z 



(39) 



The logarithmic term is present starting from the NLO {k = 1) approximation: 

while the 1/z term is present from the LO {k = 0) approximation: 

E,{as) = as Ef^ + E^^^ + a] Ef\.. . (41) 

2.5 Monte Carlo form of DGLAP equations 

The z = 1 singularity in Eq. (|3ip needs special treatment in the Monte Carlo formulation 
of the DGLAP equations. Rewritting Eqs. ^ for the parton distributions multiphed by 
X, denoted in the matrix form as 

jxgi, . . . xq^,...,xq^^, xG| (/i, x) = xD(/i, x) = Q(/i, x) , 

we have 

1 

— ^Q(/x,x)= dzP{as,z)Q(^fi,^^ , (42) 
n ^ J z 

X 

where P is the matrix of the splitting functions, which can be easily read off from Eqs. ©• 
Based on the results of Section 12. H| we can write the general structure of the splitting 
functions in the following form 

P(a„ z) = + B(«,) S{1 -z) + P{as, z) (43) 

where A, B, P are computed in powers of as- The function P{as, z) may contain singular 
terms in the limit z — > 1, proportional to powers of ln(l — z). 

For simplicity of the notation we suppress the //-dependence of the parton distribution 
and the splitting functions in the following. Substituting fjl3| to Eq. (|42| and using 
definition ()32|1 . we find 

9 Q(x/z)-Q(x) , - 



din jj 



Q(x) = Jdz!^A Q(^/^^)^ ^Q(^) + p(^) Q(^/^) 



2 

X 



{A ln(l -x) + B}Q(x) . (44) 



Now, we introduce a small cutoff in the upper limit of the integration, 1 —>■ (1 — e), which 
isolates the z = 1 singularity. Performing the integration, 

l-e 

j dzA^^^ = {Alne - A ln(l - x)} Q(x) , 



we find for Eq. 



l-e 1 



^ Q(x) = jdzA^^^ + j dzP{z)Q{x/z) 



din fx 



+ {Alne + B}Q(x). (45) 

Inserting back tlie yU-dependence, tlie equation above can be written as 

1 

a 



din fT 

witli tlie kernel 



Q(...)^/d.P(......)Q(,../.) (46) 



P(«„ z, e) = Q{l-z-e) + {A(«,) Ine + B(a,)} 5(1 - z) + P(«„ z) . (47) 

This form of the DGLAP equations is a starting point for the Monte Carlo generation. 
Let us notice that the presented formulae are valid for both representations of the parton 
distributions and splitting functions discussed in Section 12.21 Explicit expressions for the 
splitting functions up to the NLO are given in Appendix 1X1 

3 Markovian algorithm for parton distributions 

In the following we show how to transform the QCD evolution equation (DGLAP type) 
into an integral homogeneous equation and solve it by means of iteration. The general 
properties of the evolution equations and the related diffusion equations are discussed in 
Appendix O using simple environment of the discrete space. Below we discuss a more 
complicated case of the mixed, discrete-continuous, space of the QCD evolution equations. 

3.1 Classic iterative solution 

Introducing the variable 

t = ln/i = lng, (48) 
the evolution equations ()42|1 in the component form read 

d 



DK{t,x) = J2i'^Kj^Dj){t,x) 



dt 

J 



I ^7Kj{t,z)Dj(t, 

X 



(49) 
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Notice that due to the definition of the evolution variable t, the splitting functions Tkj 
are related to those from section [TJ Eqs. (jH), by 



?Kj{t,z) = 2PKj{a,{t),z) (50) 

with a possible dependence of as(t) also on z in order to accommodate coherence effects. 

The next important step is to introduce the infra-red (IR) cut e and to isolate a part 
of the kernel T diagonal both in the parton (flavour) index and in the z-variable^ 



yxAt, Z) = -y^KKit. 5(1 -Z)+ ?lj{t, Z) 

ylj{t, z) = ?K.j{t. z) 9(1 - ^ - e(t)) Q{z - e'), 



where we also introduced a facultative lower limit e' on 2;-variable, equivalent to the 
minimal global x of the evolution. Note that in the DGLAP case there is no reason for 
the IR regulator e to be t-dependent. However, for applications in the context of parton- 
shower algorithms and of the CCFM equations [17] it is worthwhile to keep this option 
open. In any case we always assume that e and e(t) are small. A more detailed discussion 
of the LO kernels is given in Appendix ^ After the above splitting of the kernels the 
evolution equation becomes inhomogeneous 

^^^(t, x) + y^^(t) DK{t, x) = X](yL ® Dj) {t, x). (52) 



It is easily made again homogeneous 

dt 



e-*K(Mo)|. (e*-(*A>)D^(t,x)) = J2i^%j^Dj){t,x), 



t 



(53) 



$x(t,to) = J dt' T^^^(t',e(t')), 

to 

and turned into an integral equation 

t 

e'''<^''''^DK{t,x)=DK{to,x) + J citie*-(*-*")5^(T|^®Dj)(ti,a;). (54) 



to 



Its another equivalent form, which is more convenient for iteration, reads 

t 

D^(t,x) = e-*-(*'*«)D^(to,x) + j rftie-*-(*'*^)5^(T|^®Dj)(ti,x). (55) 



to 



^AU components proportional to (5(1 — z) reside in the diagonal part of the matrix y^j anyway. 
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The iteration of the above equation provides a solution in terms of a series of integrals 



oo n p „ „ 



ri=l Kq,...,K„-x 1=1 



to 



X e 



i=l 



i=l 



(56) 



where fc„ = k. At this point we have many options for the MC implementation of the 
multidimensional integrals given by the above expression. Quite generally, they can be 
divided into Markovian and non-Markovian groups of the MC implementations. In the 
following we shall describe solutions of the Markovian type. However, it will be done such 
that the mechanism to switch to a non-Markovian method will be as easy as possible. 



3.2 Markovianization 

Contrary to the evolution of the non-singlet PDF or of the singlet one-component PDF, in 
the most general case represented by the Eq. (j56|) one cannot express its integrand as an 
exact product of the Markovian single-step probabilities, each normalized to 1. However, 
the general iterative solution of the evolution equation in Eq. (|^ can be expressed in 
terms of the (unnormalized) transition density 

n{ti,Xi,Ki\ti_uX,_i,Ki_^) = e{ti-t,_i) ?%^j,^^^{U,Xi/xi_,) e-*^'.-i (57) 

as follows 



oo „ n r- p p 

DK(t,a;) = e-*^(*'*«)Dx(to,x) + 5^ J2 dxo H dU dz, 

n=l Ko,...,K„.ii i=l L/ i 



to 



-'S>Kit,tn) 



(5^ 



X 



JJ n{ti, Xi, Ki\ti_i, Xi_i, Ki_i)5 {x-XoYl Zi) DKoito, Xq). 



i=l 



1=1 



The above expression looks almost as a product of Markovian transition probabilities, 
except that Q lacks a proper normalization 



oo 1 



j dti j dziY,^ituXi^Ki\U_i,Xi.i,K,^i) = y"d(Ti^,_,(ti,ti_i)) e-*^-i V 1, 

n n 



ti-i 

where 



1 



TK{tM) = j dt' j dz 5^yf^(t', 



z . 



(59) 
(60) 



to 
11 



The above problem cannot be cured by changing integration variables or normalization of 
the PDFs. On the other hand, one can define a properly normalized transition probability 



uj{ti, Xi, Ki\ti_i, Xi-i, Ki^i) = 6(ti-ti_i) '?fc^K,.^iti,Xi/xi-i) e ^\ 

oo 1 



(61) 



U-1 

and express 

Q{ti,Xi,Ki\ti_i,Xi_i,Ki_i) = e^'^^-^'-*"*'-^^u{ti,Xi, Ki\ti_i,Xi_i, Ki_i) , (62) 

where 

t 1 

A^(t, to) = TK{t, to) - <^K{t, to) = jdt' jdzJ2 '^Mt', z), (63) 

to 

which is independent of the IR regulator e{t). This opens a way to the Monte Carlo 
algorithm with weighted events, in which the Markovian algorithm is based on the u 
distributions and the correcting weight w = Yl{Q/uj) brings back the MC distributions to 
the original ones. 

As usual, Markovianization cannot be accomplished without adding one extra inte- 
gration variable. We do this starting from the identity 

oo 1 

j dti j dzi^u{tuXi,Ki\ti.i,Xi.i,Ki^i) = e-^''^-^^'^'^-'\ (64) 
t 

and then we add the (n + l)-th "spill-over" variables in the integrals using 

oo 1 

gA^jMn) j ^^^^^ j ^^^^^ Ujitn+l,Xn+l, Kn+l\tn, X„, K^) = e"*^" . 



(65) 



t 



Summarizing all the above discussion, we transform Eq. into a new equivalent 
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form 



x) 



ti>t 

1 



+ ^ / dxo / dtn+idZn+1 ^ / rfijd^;, 

n=l t„+i>t ^"+1 Ko,...,K„-i i=it.^t 

n 

n e^^-i(*-*'-^^cu(t,, Xi, Ki\U_,, Xi_i, Ki_,) 



(66) 



n 

X 

1=1 



S{x - XoY[zi) DKo{to,Xo)- 



X 

i=l 



In the Markovian Monte Carlo algorithm implementing exactly the above series of the 
integrals we neglect primarily the factor 

n 

w = e^^" Yl e^^-i , (67) 

i=l 

such that the whole series of integrals can be implemented readily as a Markovian chain of 
steps with the normalized transition probability u!{ti, Xi, Ki\ti-i, Ki^i) for each single 
step. The original integrals and distributions can be recovered by means of applying the 
MC correcting weight w defined above. The only technical problem is that w > 1 and 
one may struggle to find the maximum weight, in order to turn weighted events into 
unweighted ones. It is an unavoidable price to pay in this method. 

In the case of the single-component evolution (singlet or non-singlet) we recover auto- 
matically the constant-weight algorithm 

n 

^ ^ gA(M„) Y[ e^(**'*^-i) = e^(*'*°) . (68) 

i=l 

In the case of the non-singlet evolution we even have w — 1. 



3.3 Generation of a single Markovian step 

The description of the Markovian algorithm of the previous section is incomplete without 
providing at least one method to generate exactly the distribution of a single step forward, 
{to, xq, Ko) {ti, ziXo, Ki), in the primary Markovian algorithm 

duih, zixo, Ki\to, xo, Ko) = e(ti - to) y|,i^o(ti, zi) e-^'''>(t^''°Utidzi. (69) 
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The natural method of generating the above 3- dimensional distribution (including one 
discrete variable) can be read from the reorganized normalization integral 



oo i 

1= dti^ / dzi uj{ti,ziXo,Ki\to,Xo,Ko) 

^1 

1 „ 1 



J ^ ' E Jdz' nK^h, z') j / dz' n.^^it,, z') 

X=q,g 

1 1 

= j dr{t{) J2 PiKilh) j dz^pizi\K^,h), 

^1 

where the two final integrals and the parton sum are each equal to 1 separately. 

One may generate the first variable ti by inverting the cumulative distribution r{ti). 
Because of the possible t-dependence of the coupling constant and the cut-off parameters, 
this requires inverting the distribution r{ti) numerically or preparing look-up tables for 
TKitijto) form factors and their inverse, for each parton type K separately. 

Knowing ti, one can generate the parton type Ki according to the probability ttriKq 
proportional to / dz ^^^/^^(ti, -2)- Look-up tables of the ti dependent hkiKq branching 
ratios are needed for better efficiency. 

Finally, knowing ti and Ki one can generate the variable zi according to the probability 
distribution proportional to ^^l^^i^ol^i, ^i)- Here one can generate Zi starting from some 
approximate distribution and execute an internal rejection loop with the correcting weight, 
about which the external part of the MC algorithm knows nothing. 

As one can see, u{ti, ZiXq, Ki\tQ,xo, Kq) can be generated exactly. However, because 
of the need to pretabulate the form factors and the branching probabilities there will 
always be some irreducible numerical bias in the MC results. This requires some extra 
effort to control quantitatively and reduce, if necessary. 



3.4 Weighted Markovian algorithm 

The above Markovian scenario is close to what is used in the standard parton-shower MCs. 
Here we shall describe another class of MC solutions for Eq. ()56|) . We shall stay within 
the class of the Markovian algorithms, but our aim will be to use a MC implementation 
which allows for easy and quick transition to constrained Markovian algorithms. Quite 
generally, in the MC algorithm described above, there is a tendency of "micromanaging" 
the generation of the component sub-distributions (i.e. u distributions) such that they 
are generated exactly and there is only one extra global MC weight of Eq. (jU7j) . This is an 
efficient method but the efficiency comes at a price of using look-up tables for generation 
of the cu-distributions. 

The alternative (implemented in the MC program EvolFMC [11]) is to simplify intelli- 
gently the kernels, phase space boundaries, coupling constant, etc., such that all compo- 
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nent distributions in the MC algorithm are easily generated. The compensating weight is 
applied at a later stage to exactly retrieve the original distributions and integrals. This 
also comes at a price because an extra weight will lead to a wider weight distribution and 
a less efficient algorithm, especially if one wants to turn weighted events into unweighted 
ones at the end of the MC generation. These negative aspects can be minimized by a 
better choice of approximations and the use of internal rejection loops, wherever possible. 
On the positive side, there is no need to deal with annoying procedures of controlling 
quantitatively the numerical bias due to the use of the look-up tables. Moreover, since 
the approximate distributions reflect well the singularity structure of the integrand, we 
have better insight into the physics and a better chance to move away from the Markovian 
algorithm (if we find it profitable for some other reasons). 

Looking into the LO and NLO evolution kernels in QCD, one can see that they all 
have the following structure 

PlK{t,z) = SiKAKxit) + 6{1 - z)5iKBKK{t) + -CiK{t) + DiK{t,z), (71) 

[i — z)+ z 

where Djk{z) is completely regular. The coefficient constants A^k, Bkk,Cjk and the 
coefficient functions Djx{z) can be decomposed into the LO and NLO parts: 

4 _ 4(0) , /«^WV4« R _ o^sit) (0) , fas{t)Y ii) 

AKK[t) - ~1^^KK + ) ^KK^ BKK[t) - ~^BkK + (^~^ ) ^KK 

(72) 

C..(t)^^Cf) + (^) DMt,z) ^ ^DS(.) + (^)^<(.). 

Once we have made the above decomposition, we may readily express all the form factors 
and constants entering into our integrals of Eqs. and (jSHl)- 

The virtual diagonal IR-divergent elements in the kernel matrix and the corresponding 
form factor read as follows 

'^KKit) = 2 (^AKKit) In ^ - BKK{t) 

} } f I ^ ^^^^ 

$^(t,to) = J dt' ?^^(0 = j dt'2 \^AKK{t')ln-^ - BKKit') 

to to 

The integrated real-emission off-diagonal elements needed to generate the parton type in 
the Markovian MC is now 

1 1 

vr/x(t) = J yfAt, z) = 2 5iKAKK{t) 1^ ^ + In ^ + j dz DjK^t, z) . (74) 
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The real-emission form factors T^, K = q,g, are given by 



T^XK{t') 

(75) 

= I dt'2 \AKK{t') In + ^ CxK{t') \n^ + J2 [dz DxK{t\ z 

Finally, the form factors Ak, needed in the final MC weight for the correct overall nor- 
malization, read 



t t 
Ax(t,to) = J dt'J2 J dz 9xK{t')= J dt'!^-?'j,j,it') + J2^xKit')^ 

to ^ to ^ 

t 1 
= /"dt'2 5xi^(t') + 5^Cxi^(t')ln^ + 5^ / dzDxK{t',z 



(76) 



Having expressed all the elements in Eq. ()66|) of the standard Markovian algorithm, 
let us construct an alternative MC Markovian scenario starting from Eq. ()56|) (before 
Markovianization) . First, we simplify the kernel matrix elements 



9f^{t, z) ^ yf^(to, z) = e{z - e') 0(1 - ^ - e) 



^ X 4(0) I 1 ^(0) , A 



(77) 



where D is replaced by the constant D, which is chosen to be zero when A^^\ B^^^ are 
nonzero or equal the maximum (positive) value of Df'J^{z); see Appendix |B] The above 
simplification is of course compensated by the MC weight 



Wp 



n 



(7t 



i=i y®x,i^,_i(io, Zi) 

Let us remark that the replacement as{ti) otsi^o) of the running coupling which stands 
in front of the LO kernel might cause poor overall MC efficiency. This problem is addressed 
separately in the next section. 

The above reorganization leads us to the following new formula 

1 „ i 1 



n=l Ko,...,K„-i Q 1=1 ^ 



dZi 



X e 



n 

i=l 



n=l Ko,...,K„. 



'^K,K,^i{to, Zi)e 



DKo{to,Xo)S{x - XoY\_Zi)wp, 



i=l 



(79) 
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completely equivalent to Eq. (j^Bj) . Markovianization is now done for the variant of the 
above formula in which wp is neglected. We define a new transition probability as follows 



oo 1 



U-i 

where 



U 1 



fK{u,u^i) = j dt' j dz Y,y%{t\ 



as{to 



(0) 
XK 



X 



X 



ti-1 

and the probability rate of the parton transition ^ / is now constant 
1 



/ 



niK = I dz 7fj^{to,z) 







TT 



1 



1 



(80) 



(82) 



independent of t; see Appendix El for explicit formulae. Summarizing, the transition 
probability 

oo{ti,x,,K,\U_i,x,_i,Ki^i) = e(ti -t,_i) 'J'tK,_,{to,x,/x,_i) e-(*»-*-i)^^.-i (83) 

is now such a simple function that can be generated using elementary MC methods, 
without any pretabulation. The last thing necessary, as usual, for the Markovianization 
is introduction of the "spill-over" variable. This is done with the help of the identity 



-*X„(i,tn) 



gAK„(tA) / / ^^^^^ ^ u{tn+l,Xn+l,Kn+l\tn,Xn,Kn) 



K. 



n + l 



where 



Ak{U, tj-i) = TxiU, ti_i) — ^k{U-, U-i) = {U — U-i)Rk — ^k{U, U~i) ■ (85) 



Let us stress that now, contrary to the previous standard Markovian scenario, A has 
an explicit residual dependence on the IR cut e which is necessary to cancel exactly the 
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analogous dependence of the average weight wp (similarly as in a typical MC algorithm 
for QED exponentiation). 

Summarizing all the above discussion, we transform Eq. (|79|) in a new equivalent form 



ti>t 

1 



dtidzi 



+ ^ dXo / dtn+ldZn+l ^ ^ Y\ 

n=l K„+i Ko,...,K„-i i=\^t 

n 

1=1 

n 

5{x- XoYl Zi) DKoito, Xq) WpWA , 



(86) 



X LUi 



X 

i=l 

where 

n 

Wa = e^^" JJe^^-i(*' . (87) 

i=l 

In the MC generation we proceed as before. Neglecting the weight w = wpWa, we 
generate primary MC events using the Markovian algorithm with the simplified transition 
probability a). The original distributions and integrals are recovered by applying the 
correction weight w. As already stressed, the MC efficiency will be worse than in the 
previous case, but the whole MC program is now much simpler and most likely provides 
better control of the technical precision. 

Let us note that we shall still need a precise 1-dimensional pretabulation of all the 
form factors $i^(t,to), entering into wa through Ak- 



3.5 Importance sampling for running as{t) 

In the above we took into account the t-dependence {t = InQ), that is running, of the 
strong coupling constant as{t) by reweighting MC events. This is very inefficient and it is 
rather easy to introduce the relevant t-dependence of a<j(t) at least at the one- loop level 

already in the underlying MC distributions. 

One can see that in the i-integration in Eq. (f?^ we have effectively 

/27r 271 C 
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It is therefore natural to introduce the variable r^, 



Ti = ln{ti - InAo), U = In Ao + exp(ri), 
dti = {ti — In Ao)(irj = e'^^dTi, 



(90) 



instead of'^ ti. This change of variables will lead to the Jacobian factor e"^' = ti — In Aq in 
the integrand, which will cancel the unwanted factor e""^' = {t — In Aq)"^ present in as{t) 
in the MC weight . The two-loop and more complicated contributions to as{t) may still 
be added by reweighting events, without spoiling much the efficiency. 

We start again from Eq. ()56p (before Markovianization) . The kernel matrix elements 
are now simplified more "gently" with respect to Eq. ()77|) 



?%it, z) ^ T«,(t, z) = Q{z - e') 6(1 - z - e 



36 



TT 



X 



l-z) 



A 4(0) _L ^ -u n 



(91) 



where Djk = Djk and e{t) — > e. The new compensating MC weight is 



i=l -^K^K.^^K^^^ Zt) 



(92) 



The above reorganization leads us to the following new formula, completely equivalent 

to Eq. (jnni), 



1 



\^ II - n n 

n=l Ko,—,Kn^i i=l ^ 



X 



1 

/n 
dxoYl 
i=l 



(93) 



i=l 



In the following we shall use a new function 



(94) 



because it (foes not depend on Tj anymore. This is the whole point of the ti — > change 
of variables. 



■^Of course, we are aware of a possibility of introducing r as an evolution "time" in the original 
differential equation from the very beginning. We proceed this way in order to get more insight into 
various versions of the MC algorithm. 
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Omitting wp, we proceed to Markovianization following the example of the previous 
section 



oo 1 

/*,/d..E^(n,...A-<|r.-...-.,i^.-0-i. 

Ti-l ^» 



where 



(95) 



fK{Ti,Ti.^)= / dr' dzY,"^: 



JK\ 



n-i 



in - T, 



X 

and the probability rate of the parton transition K I is now a constant 

1 



(0) 
XK 



R 



K 



TtiK = / dz yf^(;2) 



/3o 



5iKA%\\n-_ + Cy^Xn-^ + DiK 



(96) 



(97) 



again independent of t. The final transition probability to be generated in each step of 
the Markovian algorithm reads 



where 



-(ti— ri_i)-Rif^_j 



(98) 



X 



1-z 



o/jc^xii: "I — ^iK + -'^-f-f 
z 



(99) 



is again a simple function which can be generated, without any pretabulation, using ele- 
mentary MC methods. The overall recipe, as compared with the previous MC algorithm, 
is to replace: U — ^ and as{to)/T: —>■ 2/(3o in the generation of the primary MC distribu- 
tion cu, before applying ivp. 

Inevitably, to complete the Markovianization, the integral over a "spill-over" variable 
Tn+i is added with the usual identity 



00 



g-*K„(r,r.) ^ gAK„(r,r„) J ^^^^^ J ^^^^^ a;(T„+i, X^+l, K^+l^, X„) , (100) 

Kn+l 



'n+l 

T 
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where 



(101) 



The final formula, equivalent to original Eq. (f?^ . for this MC scenario with the im- 
portance sampling for the running ag reads as follows 

Di^(r,x) = e^^("'"°) j dridziY,^in,xi,Ki\ro,x,K) Dk{to,x) 

ri>T ^1 

1 

+ ^ dxo / dTn+ldZn+l JJ / dTidZi 

n 
i=l 

n 

X 5(2; - Xo J]^ ^i) Dko {to, Xo) WpWA , 



n=l 



(102) 



where 



= e 



AK„(T,r„) TT Ajf^_^(Ti,ri_i) 



n 



(103) 



The above formula looks almost identical to Eq. (|86|) . 



4 Markovian MC for parton-momentum distributions 

The factor 1/z in the bremsstrahlung kernels causes a significant loss of MC efficiency 
due to exp(Aii-) which contains uncompensated In(e'). We can get rid of this annoying 
phenomenon by switching to the xD{x) which evolve with the kernels zP{z). The reason 
for improvement is that kernels zP{z) fulfill the momentum-conservation sum rules. The 
evolution equations for xD{x) read 



i 

dtxDK{t,x) = 22 / —z?Kj{t,z) -Djit,-]. 

J z z \ z / 



(104) 



The iterative solution can be obtained from the above formulae, or equivalently by 
multiplying both sides of Eq. (fSBj) by x, 

xD,,(t,x) = e-*^(*'*o)xDi^(to,x) + ^ / c/xo n dU e{ti - U.i) dzi 

n=l { Ko,...,K„^ii=l { 



1=1 1- 



X e-*^(*'*")n ^^'^tK,-.it^^z^)e~'''''-^^'"''-'^ xoDK,{to,xo)6{x-xol[z,) 



i=l 



(105) 
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where K = Kn- It was essential to exploit the condition W^^i Zi imposed by the 

overall 5-functions''. We also feel free to introduce such an overall factor x in the xD{x), 
because our ultimate aim is to use a constrained Markovian algorithm, hence such a factor 
will be dealt in the MC separately and independently with other dedicated MC methods. 
At the technical level, we may multiply both sides of the above equation by 1/x and 
absorb 1/x in the MC weight, pretending that we generate D{x) distribution as before; 
in such a case the fluctuations of the weight in the histograms of x will change but the 
distribution of x will be the same. The main change will be in the probability distribution 
cu for the forward leap in the Markovian random walk. 

Before we enter into details of the Markovian MC, let us introduce the evolution 
variable r, similarly as in the previous section 

t 

1 f , , , dt acft/i) . s 

T = — — / dh as{h), — = 106 

as{tA) J dr as{t) 

tA 

In the above transformation we may use various choices of tA and of as{t). For instance, 
we may employ the same as{t) as in the evolution equations (LO or NLO) or we may 
stay with the one-loop LO: a^^\t) — 27r/(/?o(i — InAo)). In the latter case, with Ia chosen 
such that a^s'\tA) = 27r//9o (s-g- ^ In Aq = 1, tyi = lii(eAo)), we recover the definition 
T = \n(t — InAo) of the previous section. Let us adjust tA = to to the starting point of 
the evolution and use a^s'\t) in the definition of r. With such a choice we have 



^ 1 „ r 1 

oo 



xDK{r, x) = e-^^^^'^^^xD^lro, x) + J2 dxo H / "^^^ " ^-i) / 

n=l { Ko,...,Kn-i i=l { 

n r- -1 n 

xe-*-(-'-")J] Vlj,^_^{n,Zi)e-^-i-^^^^'^^-'^ xoDKo{ro,xo)S{x-xol[zi), 



dzi 



i=l 



(107) 



where K = K^, and 



(0)/. \ 

VIkU^, z,) = z.n.K^.M^ (108) 



which depends on very weakly. In the LO case it is completely independent of r^. 

In the general (NLO) case we may decompose the evolution kernels multiplied by z as 
follows 

\z'yiK{t, Z) = zPiK{t, Z) - ^ 5iKAKK{t) + S{1 - z)5iKBKK{t) + FiK{t, z), (109) 

^ t ~ Z) + 



*This way we could introduce any power of x, say x", in front of D(x) and z°' in the kernels. 
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where the finite part F can be expressed in terms of the previously defined functions A, C 
and D as follows 



FlK{t, Z) = zDiK{t, Z) + CiK{t) - SiKAKK{t). 



(110) 



In the MC we replace, as before, the full kernels for the real emission with the LO 
approximation with the constant IR regulator e < e{t): 



VfAr, z) Vf^ro, z) = Q{1 - z - e)^^zPf^{z), 



[111) 



The approximate kernels do not depend on r. The corresponding compensating weight is 
now ^ ^ 



i=l ' KiKi_il^O, Z^) 

The probability of the forward Markovian leap reads now as follows 

oo 1 

J dn J dzi^Lj{Ti,Xi,Ki\Ti-i,Xi-i,Ki_i) = 1, Zi=Xi/xi-i, 



(112) 



n-i 

where the new real-emission form factor is defined as follows 

1 



n 1 

J dr' J dz J2^fKiro,z) 



(113) 



= {n - Ti-i) 



TT 



1 



= {n - Tj_i) ^ TVjK = in - Ti_i) Rk- 

J 

The rate of the parton transition X — > 7 is now 
1 



J dzVfK{To,z) 



TT 



z)dz 



On the other hand, the exact virtual form factor is 



T 

^K{r,To) = J dr' 



To 



, c^f (to) 



AKK{r')\n—-^BKK(r') 



(114) 



(115) 



(116) 
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In the LO, for the one-loop af'' and if, in addition, we choose e(T) — e — const, then the 
virtual form factor becomes simply 



TT 



^K{r,To) = {T-ro)^^^[A^;^>^ln--BP^) . (117) 



Summarizing, the final transition probability to be generated in each step of the Marko- 
vian algorithm reads 

(I;(t,,x,,X,|t,_i,x,_i,X,_i) = e(T, -T,_i) P^,K,.,{ro,Xi/xi^^) e-(^^-^^-i)^-i-i , (118) 
where 

VfAro, z) = Q{l-z- ^)^^[-:^/ikA'Sk + ■ (HQ) 

To complete the Markovianization, the integral over the "spill-over" variable r„+i is 
added with the help of the usual identity 

oo 1 

g-*K.(r,r„) ^ gA;,„(r,r„) j ^^^^^ J ^^^^^ ^ a;(r„+i, Xn+i|T„, X^, K^) , (120) 

T ^"+1 

where Zn+i = and 

A;^(Tj, Ti_i) = 7V(Ti, Ti_i) - (Ti, Ti_i) = (Ti - Ti_i)^;f - (t^, Ti_i). (121) 

The advantage of the method outlined in this section is that at the LO level we obtain 
for e = e 

A;^ = (122) 

due to the fact that the kernels obey the momentum sum rules. In most renormalization 
schemes (e.g. MS) this will be also valid at the NLO level^. 

The final formula for this MC scenario with the importance sampling for the running 
as reads 

xDk{t,x) — e^'^^'^''^'^ j drid^i ^^a;(Ti, Zix, XiItq, X, X) xDk{tq,x) 



T1>T ^1 



1 



+ ^ dXo / dTn+ldZn+l ^ ^ JJ / d^dzi 

r„+i>T -^"+1 Ko,...,K„-i i=^ri<T (123) 

n 

i^.-l) 



i=l 



S (x - Xo Y\ Zi) XoDko{to, Xq) WpWA- 



X 

1=1 



small non-zero value of Afe may bo present for technical reasons, that is, if we use shghtly simphfied 
kernels at the low MC generation level. 
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where Zi = K = Kn and 



n 



_j(Ti,Ti_l) 



(124) 



5 Numerical tests 

We have performed comparisons of the MC solution of the DGLAP evolution equations 
implemented the program EvolFMC [11] with another solution provided by the non-MC 
program QCDnuml6 [3]. In both cases we have evolved singlet PDF for gluons and three 
doublets of massless quarks from Qq = 1 GeV to Q = 10, 100, 1000 GeV. The comparisons 
have been done both for the LO and the NLO evolution, including the running in the 
corresponding approximation. 

In our test, we have used the following parameterization of the starting parton distri- 
butions in the proton at Qo = 1 GeV: 

xDcix) = 1.9083594473 ■ x-^-^{l - xf-^, 
xDq{x) = 0.5 • xDsci,{x) + xD2uix), 
xDq{x) = 0.5 ■ xDseaix) + xDd{x 
xD^^^{x) = 0.6733449216 ■ x-^-^{l - x 



0.2/1 ^\7.o (125) 



0.5n ^\3.0 
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xD2u{,x) = 2.1875000000 ■ x" °(l - x 
xDd{x) = 1.2304687500 ■ x°-^(l - xf'^, 

The first results of these comparisons for the LO evolution were presented in Ref. [6] . 
They showed a 0.2% discrepancy for the gluon distributions between EvolFMC and QCDnuml6. 
This numerical bias is eliminated in this paper. In Figs. we show the resulting gluon 
and quark distributions evolved to Q = 10, 100, 1000 GeV in the LO approximation. As 
one can see, these two calculations agree to within 0.1% for the gluon as well as for the 
quark-singlet distributions. The origin of the previous 0.2% discrepancy for gluon was 
identified as a result of too high values of the dummy IR cut-offs: e = e = 10"'^. The new 
results have been obtained for e = e = 10"'*. In the small-x region [x < 0.1), we have 
found a similar agreement with the program APCheb33 [4], which uses the Chebyshev- 
polynomial technique to solve the DGLAP equations. 

In Figs. OlE] we present the results of similar comparisons for the NLO evolution. In 
the NLO, there is some ambiguity in calculation of the running a^. In EvolFMC we use 
a definition of given in Appendix^ However, in the original version of QCDnuml6 a 
different definition of at the NLO was employed. For the sake of our comparisons, we 
have replaced in the QCDnuml6 code the routine for as evaluation with the appropriate 
routine from EvolFMC. We have checked for several values of Q^, that after replacement 
the two programs give numerically the same values of OsiQ"^)- The NLO results for the 
gluon and quark-singlet distributions from EvolFMC and QCDnuml6 agree within ~ 0.1%, 
as in the LO case. Again, for a; < 0.1 we have found a similar agreement with APCheb33. 
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Figure 1: The upper plot shows the gluon distribution xDG{x,Qi) evolved from Qq = 1 GeV 
(black) to Qi = 10 (red), 100 (green) and 1000 (blue) GeV, obtained in the LO approximation 
from EvolFMC (solid lines) and QCDnumlG (dashed lines), while the lower plot shows their ratio. 

The MC calculation for the NLO evolution is, of course, slower and less efficient than 
for the LO one but not very much. In comparison with the LO results given above, the 
NLO results were obtained for the statistics about 2 times higher and their computation 
required about 3 times more CPU time. 

All the above results were obtained in the weighted-event mode of EvolFMC. In the LO 
case the weighted events can be turned into the unweighted (weight = 1) ones without 
difficulty - the event weight is well-behaved, non-negative and bounded from above. At 
the NLO, however, the situation is problematic. As was described in Section l273t the Ppc 
and Pgf splitting functions at the NLO acquire logarithmic singularities at 2; = 1. In 
our MC algorithm this leads to large positive weights for the F ^ G transitions and to 
negative weights for the G ^ F transitions in the region of 2; > 0.95. While the problem 
of large positive weights can be cured with appropriate importance sampling, there is 
no technical method to turn negative- weight events into unweighted events. The current 
version of EvolFMC implements only the weighted-event solution for the NLO DGLAP 
evolution. These problems will be addressed in our future works. 
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Figure 2: The upper plot shows the singlet quark distribution Dq) evolved from Qq = \ GeV 
(black) to Qi = 10 (red), 100 (green) and 1000 (blue) GeV, obtained in the LO approximation 
from EvolFMC (solid lines) and QCDnumlG (dashed lines), while the lower plot shows their ratio. 

6 Summary and outlook 

In this paper we have presented in detail two Markovian MC algorithms for solving the 
DGLAP equations up to the NLO accuracy. The one of them is based directly on the 
evolution of the PDFs and the other one instead of pure PDFs uses the parton-momentum 
distributions. The latter algorithm is more efficient due to the momentum-conservation 
sum rules. The evolution is done in parton flavour space, in the current version for gluons 
and three light quarks but it can be easily extended to include more flavours. Both 
the above algorithms have been implemented in the MC generator EvolFMC (written in 
C-I-+). This program has been cross-checked against two independent non-MC programs: 
QCDnuml6 and APCheb33. The numerical tests show that with today computer CPU power 
the Markovian MC is able to solve efficiently and precisely (to the per-mil level) the QCD 
evolution equation up to the NLO. Therefore, it can be used to cross-check other, non-MC 
methods or even as an alternative to them. As was pointed out in Introduction, the MC 
method is not competitive with other techniques in terms of the CPU time but is has 
certain advantages that may be important in some cases - it is usually less biased and 
more stable numerically as well as more flexible for possible extensions. 

So far we have included only light (massless) quarks in our MC algorithm, however. 
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Figure 3: The upper plot shows the gluon distribution xDa{x,Qi) evolved from Qq = 1 GeV 
(black) to Qi = 10 (red), 100 (green) and 1000 (blue) GeV, obtained in the NLO approximation 
from EvolFMC (solid lines) and QCDnumlG (dashed lines), while the lower plot shows their ratio. 

extending it to heavy quarks does not pose any problem. It can be realized either by simple 
rejection of extra massless quarks below mass thresholds or by some importance sampling 
that accounts for these thresholds. Also extension to the NNLO seems straightforward - 
one only needs to implement the NNLO evolution kernels [15,16]. A more serious problem 
concerns the divergences of the NLO kernels at 2; — > 1 which for some transitions lead to 
negative weights. This indicates that some resummation in this region might be necessary. 

As was mentioned in Introduction, we do not consider the Markovian MC algorithm 
described in this paper to be only a tool for solving numerically the evolution equations for 
parton densities. It can be also used as the basis for constructing the FSR parton shower 
MC event generator which would generate physical events in terms of particle flavours 
and four-momenta. This algorithm cannot be directly used for the ISR evolution (parton 
shower), since in that case one needs to impose some energy- momentum constraints on 
the partons entering the hard process. However, it can be used as a starting point for 
developing any kind of constraint MC algorithms as well as it can play a role of a testing 
tool for corresponding MC programs, see Refs. [7-10]. 

Yet another way of development of the above Markovian MC algorithm can go into 
the direction of solving the CCFM equations [17]. In this case one deals with the so- 
called unintegrated parton distribution functions, which in addition to x and depend 
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Figure 4: The upper plot shows the singlet quark distribution xDq evolved from Qq = 1 GeV 
(black) to Qi = 10 (red), 100 (green) and 1000 (blue) GeV, obtained in the NLO approximation 
from EvolFMC (solid lines) and QCDnumlG (dashed lines), while the lower plot shows their ratio. 

on the partonic transverse momenta kx- This may be important for better modehng 
physical events in deep inelastic scattering as well as hadron-hadron collisions, see e.g. 
Ref. [18]. We have already implemented the so-called one-loop approximation of the 
CCFM equation [19-21] which will be reported in our forthcoming paper [22]. 

A QCD kernels at NLO 

A general parton-parton transition matrix for a gluon and three quark flavours {d, u, s) 
can be written as 





-G, 


Pg. 


-d, 


Pg^u, 


Pg^s, 


Pg. 


-J) 


Pg^u, 


Pg. 


—s 


Pd. 


-G, 


Pd. 


-d, 


Pd*—ui 


Pd^s, 


Pd. 


-J) 


Pd.—ui 


Pd. 


-s 


Pu. 


-G, 


Pu. 


-d, 


P 


P 


Pu. 


-J) 


P - 


Pu. 


-s 


Ps. 


-G, 


Ps. 


-d, 


p 


p 


Ps. 


-J) 


p - 


Ps. 


-s 


Pd. 


-Gi 


Pd. 


-di 


P- 


Pd^si 


Pd. 


-J) 


P- 

^ d^m 


Pd. 


-s 


Pu. 


-G, 


Pu. 


-d, 


P- 

n<— n; 


P- 

n<— S) 


Pu. 


-J) 


P- - 


Pu. 


-s 


P-s. 


-G, 


P-s. 


-d, 


P- 


P- 


P-s. 




P- - 


P,. 


-s 



(126) 
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where Pj^i = Pj^-i{as, z). At the NLO, the kernels can be decomposed as follows 



P{as,z) 



27r 



P(°)(^) + 



27r 



PW(^), 



(127) 



where the NLO QCD coupling in the M5'-scheme is 



a,{t) = af\t) |l - af\t) ^ ln(2[t - InA^])} , 
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and t = InQ (n/ is the number of active flavours). 
The LO kernel matrix takes a simple form 



p(0) 


p(0) p(0) p{0) 
GF; GF; GF; 


p(o) p{o) p(o) -1 

GF; GF; GF 


p(0) 

p(0) 
^FG' 
p(0) 
^FG-i 


p(0) Q Q 
FF' 

P^°^ 

0, 0, pfl 


0, 0, 
0, 0, 
0, 0, 


p(0) 

^FG^ 

p(0) 

^FGi 

p(0) 

^FGi 


0, 0, 0, 
0, 0, 0, 
0, 0, 0, 


P^FF^ 0, 

0, pS, 
0, 0, pfl _ 



where 



p(0) 
^GG 



[z) = 2Ca 



(1^-2 + ^^^-^) + .. 



+ S -8{l-z), 



f ) (z) = [z' + (1 - z)^ 
1 + (1-^)2 



gf(^) - 



(0) 
FF 



(z) = C, 



z 

[ 1 + z^ 



(128) 



(129) 



(130) 



and the colour-group factors are: Ca — Nc — 3, Cp 
The NLO contribution to the kernel matrix can 



= {Nl - l)/2N, = 4/3, Tr = 1/2. 
be expressed in the following form 



P«(^) 



p(l) 

^GGi 


p(l) p(l) p(l) 
GF' GF' GF; 


p(l) p(l) p(l) 1 
GF; GF' GF 


p(l) 

^FGi 

p(l) 

^FGi 

p(l) 

-^FG^ 


pV/+6'(l) p6'(l) p5(l) 

p5(l) py+S(l) p5(l) 

p5(l) p5(l) py+5(l) 
59 , -t gq , J^qq , 


pV/+A'(l) pi'(l) pi'(l) 
qq ' 99 ' 99 
p5(l) pV+.^Cl) p5(l) 

99 ' 99 ' 99 
pS(l) p5(l) pV+S{l) 
^qq 1 ^qq i ^qq 


p(l) 

^FGi 

p(l) 

^FG^ 

p(l) 

^FGi 


pV/+6'(l) p6'(l) pS il) 

qq 1 qq ' qq ' 
5(1) pV+S(l) p5(l) 
qq 1 qq ' 99 ' 

pS(i) p5(i) py+s(i) 

99 ) 99 ) 99 ) 


pV+S(l) pS(l) p6'(l) 
-Tqq , rqq , rqq 

pS{l) pV+S{l) pS{l) 
rqq , rqq , rqq 

pS(l) pS(l) pV+S(l) 
rqq , r-qq , r^qq 



(131) 
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where P: 



(1) 



p]jj{z) and we have used a short-hand notation P^j 



ij 



V+S{1) 



pV{l) p5(l) 



IJ 



The above matrix can be simphfied by exploiting the identity 



p5{l) ^ pS{l) 
^qq ^qq i 



(132) 



which is true up to the NLO. However, we prefer to keep a more general form of the 
kernel matrix which can be useful for some tests and possible future extensions. The 
above kernel matrices can be easily extended to include more quark flavours. 

The non-singlet and singlet-quark kernels can be expressed in terms of the basic NLO 
splitting functions P+, P_ and Pff-, defined in Refs. [13,14], as follows: 



pV{l) ^ ^ 



p(l) + p(l) 



pV(l) 



(1) 



P 



(1) 



1 



2n 



pPf - 



. (133) 



Finally, all the elements of the above kernel matrix can be calculated (e.g. numerically) 
from the six basic NLO splitting functions of Refs. [13, 14] 



p(i) p(i) p(i) p{i) p(i) p{i) 



(134) 



These basic splitting functions are given at the NLO by the following expressions: 



Pi'\z,e) 



Q{l-z-e) + bI\z) + 



C 



(1) 



Af Ine 



5{l-z), 



(135) 



P'^l{z,e) = -^Q{l-z-e)+B','\z) + 



A 



(1) 



A 



(1) 

G 



1 - 2 



e(l-^-e) +Pg)(z) 



C« + 4hne' 



C^+Aghne" 



5(1 -z), 
5(1 -z). 



(136) 
(137) 



p''pI{z) = -CfTr<| 4-9^+ (4z- l)ln2 + (2^- l)ln^z + 41n(l - z) 



+ 



+ -CaTr 



10 - -TT^ + 21n2 

182 14 40 

\ z H 

9 9 9^ 



1 -2 



Z 



-4 In 



1 - z 



z 



+ (1 - ^Y] } 



136 38,, ^, ,^ , 
z ^ I m 2; — 4 m(l — z) 



3 3 

(2 + 8z) In^ z + 2 [z^ + (1 + ;2)2] ^2(^) + 



'1381 



71^ 218 

y ~ ~9~ 



44 

-I In 2: - In^ 2: + 4 ln(l - z) - 

3 



2\n\l-z)\ [z' + il-zr]^ 
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P^%) -Cf[-1-1^- 2^1n(l -z)- 1±1^^ l„(i _ ^) [3 + 1„(1 - z)] 



where 



+ ( 2 + 1 In^ - ( 1 - 1 In^^ 



1 + (1 - z)^ 



+ 2^1n(l - z) + 



1 11 



+ ln2(l - z) -2\nz ln(l - + -In^;^ 

2 



l + + 



CpTfi-z + 



1 + il-z) 



^5 

^(1) _ ^2 



2 

67 _ TT^ 



1 

2 
20 



17 llvr^ 



60 



C^Tf ( . 



1 27r2 
6 ^ ~9~ 



9 



Tf 



-Ca + Cf Tf, 



67_ TT^ 

9 3 



20 



with Cs = C(3) ~ 1.2020569. 

The non-singlet coefficients take the form 

A^±\z) = Cl [-2{l+z^)\Yiz\n{l -z)-Zhiz\ 



+ -CfCa 



11 



67 7t' 



In ^-|- — ln2;-|-— 



9 



(1 + -^) 



In^; + 



(1 + ^'), 



B^\z)^Cl 



-2zlnz--{l + z)ln'^z-5{l-z) ± Pa{z) 



40 

2(l + z)ln^ + -(l-z) T Pa(z) 



C^Tf-il-z), 



where 



Pa{z) = 2 S2{z) + 2(1 + ^) In ^ + 4(1 - z) . 

1 + z 
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Notice that, as it should be, A^.\l) = A^g\ The function 5*2 is defined in the interval 
(0,1] 



1/(1+^) 



S2(z)^ I ^In^^ - -2Li2(-^) + ^In^z - 21n^ln(l + ^) - ^, 

J y y 2 6 



z/{l+z) 

where the dilogarithm function Li2 is given by 



(144) 



(145) 



with a branch point discontinuity on a complex plane along [1, +oo]. 
For the singlet and gluon coefficients we have 



1 + ^ + - [1 - 3^ - (1 + 2;) In 2;] In z 



- 2 



14^ 
\-z 



+ ln(l - z) 



ln^+ 2 



1 



1 + ^ 



+ Cf Ca 



67 7r2\ ^ 14,^ 

(1 +^) H (1 - z) 

18 6 / ^ ^ 3 ^ ' 



+ 



1 fl + z' 



2\l-z 



— + lnz \nz- — — 

3 / 1 + z 



^ ^ ( 16 10,^ , 40 



S2iz) 

40 40z 112^2 2 fl + z^ 



9 



3 V l-;2 



In 2; 



+ ^2 + 10 z + y ^2 ) In z - 2 (1 + In^ 2 



16 



(146) 



+ 4(1 + ^) In^ ;2 + [In^z - 41n2;ln(l - z)] 

-L ^ 



+ 



In^ z - 4 In 2; ln(l - z) + — - — 

9 3 



^-2 + z{l-z) 



+ 2 ( - - - 2 - z - ) ,52(^) 

1 + 2; 2; 



+ C^r^^2(l-^) + ^(^2-^-^)-l(l + ^)ln^-^ 



-2 + z{l-z) 



r 20 4 1 

+ Cj. T J 8(^ - 2) + —z^ + _ - (6 + 10^) In ^ - 2(1 + ^) In^ ^ L 



(147) 
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It is easy to check that the above functions have no singularity at z = 1. 

It is worth noticing that the non-singlet quark kernels have particularly simple ana- 
lytical representations 



Alpine 



S{l-z) + 



e{l-z-e) 



-2z\nz- -(l + z)!!!"^ z-5(l- z) 



(148) 



+ -CfCa 



40 

2(l + ^)ln^ + — (1-^) 

o 



-Cj.Tf-{l-z), 



P^^'\z)^Cf{Cp-^Ca]Pa{z). 



(149) 



Using these analytical formulae in numerical evaluations of the kernels P^q^^^ and P^ 



>v(i) 



can be faster and more stable numerically than computing them indirectly from the split- 
ting functions 

In the Monte Carlo implementation, the above kernel matrices are used for generation 
of real-parton radiation, i.e. for z < 1 — e. This has to be compensated by the appro- 
priate Sudakov form factor summing up virtual and soft-parton corrections, i.e. terms 
proportional to 5(1 — 2;). At the NLO, the Sudakov exponent $x takes the form 



t2 

^K{t2,tl)^2 J dt 



O^sjt) 5(0) 

27r 



2n 



(150) 



tl 



where the NLO as{t) has to be taken for both terms. The factor of 2 in front of the integral 
is due to the fact that our evolution "time" ist — \n.Q. Integrating over r = ln(i — In A), 
we obtain for gluons 



$G(r2,ri) 



/So 



[Xo(r2) - Xo(ri)] 
+ [Xi(t2) -Xi(ti)] 



2CA\Tl--^CA + \Tf 

e D 6 



(1) 



(151) 



and for fermions (both quarks and anti-quarks) 

2 



$f(t2,Ti) 



/?0 



2CfI^-- -Cf 
e 2 ' 



where 



[Xo(t2) - Xo(ti) 
+ [Xi(t"2) -xi(n) 



(152) 



(153) 
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and 



r2 + 2 ln2 + 



r + - I + In^ 2 



-r-ln2 - - 1 e~" + 1 



(154) 



B Gluon and quark-singlet kernels at LO 

As already advocated, we isolate all singular parts from the kernels (cf. Eq. ((71])) 

PlK{t, Z) = —^-—6iKAKK{t) + 6{1 - z)5iKBKK{t) + -C/^W + DiKit, z). (155) 

(1-2;)+ z 

After expanding to LO (see Eq. ()72|) ) the resulting components A^'^\ B'^^\C^'^\ D^'^'> are 
listed in Table [T] 



IK 


.(0) 


r(0) 


MO) 




Dik{z) 


JdzDf^iz) 


GG 


ICa 


11 /~i '2rp 


2Ca 


2Ca{-2 + z-z^) 
2Tfiz' + il-zr) 
Cf{-1 - z) 







qG 
qq 


2Gf 








2Tf 



ITf 


Gq 






2Cf 


Gf{-2 + z) 





-|Cf 



Table 1: The elements of the singlet LO kernels {Tf = njTji). 



Let us consider evolution of the simple two-component state consisting of the gluon 
and the (singlet) quark with the LO evolution kernel 



p{0), 



pSkz), pS^z) 



Pc 



where 



z) 



(156) 



(0) 



FF ' 

3(0) 
GF 1 



(157) 



and Pqqi PpGi PpFi Pgf given explicitly in Appendix 1X1 

Let us concentrate now on the second (less standard) Markovian algorithm described 
in Section 13.41 The definition of the simplified kernel matrix elements, Eq. fl77|) 



'yfj,{to,z) = Q{z-e')Q{l-z-e 



n 



1 - z 



A 4(0) , ^r^W , n 



(15^ 
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needs Dik^ which are also provided in Table [TJ 

In this simple case we may explicitly list the parton K I transition rates 



TT 



(159) 



as follows 



TT 



2Ca[ Ini + lni , 2C^lni 



2T, 



2Cf In 4 



(160) 



and also the characteristic decay rates Rk = Xlx ^xki K = g,q {in a. primary MC) 

Ra='^i2cJln'^ + ln^]+2T, 



R„ 



TT 
TT 



|2C^ln^ + 2Ci.ln^|. 



(161) 



C LO kernels for parton-momentum distributions 



IK 


AO) 






maxF)^(z) 


jFi'^{z)dz 


G^G 
q^G 
q^G 


2Ga 




11 2rp 





2Gaz{-2 + z-z^) 
Trz{z' + (1 - zf) 
Trz{z^ + (1 - zf) 




Tr 
Tr 




G^q 

q^q 
q^q 




2Gf 





|Cf 




Gf{2 -2z + z^) 
Gf{-2 -z-z^) 



2Gf 







G^q 
q^q 
q^q 





2Gf 





^Gf 


Gf{2-2z + z'^) 


Gf{-2 -z-z^) 


2Gf 




^Cf 


17 r' 



Table 2: The elements of the LO kernels for parton-momentum distributions (Ty = iifTR). 

We may decompose the evolution kernels for parton-momentum distributions in the 
LO and NLO as follows, see Eq. dTUnD 



zPiK^t, Z) 



l-z) 



-5lKAKK{t) + z)5iKBKK{t) + FiK{t,z) 



(162) 



where q represents one quark flavour and the corresponding LO components are listed in 
Table |21 Comparing to the previous appendix, we have 



FlK{t, z) = zDiK{t, z) + GiK{t) - SlKAxKit) 



(163) 
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In the MC we use the simphfied kernels (Eq. 

,(0) 



VfAro, z)=eil-z- ^)'^{y^JikAP, + (164) 
In this simple case we may explicitly list the parton transition K ^ I rates (Eq. (jlisp ) 



1 



7r/x(ro) = / dzVfKiro,z) = 



TT 



e 



(165) 







where 



1 

ffi^ jdzFf^z) (166) 




are listed in the Table |21 

The momentum-conservation sum rule reads 



^ j dz zPxxit, z) = BKK{t) + ^jdz FxK 

^0 ^0 



(t,z) = 0, (167) 



which is manifest in Table |2] for the LO case but it also holds for the MS NLO kernels. 
Note that the above sum rule determines unambiguously the virtual part of the kernels 



1 



BKK = -Y,f^J< = -Y. Pxk{z), (168) 







X X 

both in the case of the LO and the NLO. 

D Generic discrete Markovian process 
D.l Diffusion and evolution equations 

Let us consider a general "evolution equation" for the multistate discrete system 

dtNi{t) = Y,PiL{t)NL{t). (169) 

L 

Let us ask whether the time dependence of the above system can always be interpreted 
(implemented) as a probabilistic stochastic Markovian process, i.e. in terms of MC events 
with weight equal to 1. We shall see that this is not true in the general case and we 
shall show under which restriction on the transition matrix Pn the above conjecture on 
Markovianization is true. Weighted MC events are excluded from the consideration, i.e. 
by the Markovian process we understand the probabilistic process with weight=l events. 
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Without any loss of generality, in our starting Eq. ()169p we have chosen a discrete system 
in order to simplify the reasoning. 

An answer is found by examining a general "diffusion" process in the discrete space. 
We shall derive the evolution equation ()169p finding restriction on the transition matrix 
PiL- Following this path of reasoning we first define a general transition probability of 
an object which is exactly in the state / at the initial time to- The explicit transition 
probability {to, I) — >■ {ti,K) at the later time ti > to into another state K ^ I is defined 
as follows 

-fdt'ZPjiit') 

dp{K,ti\I,to)\^^, = 9{ti - to) Piuih) dti e '0 -'^^ . (170) 
It is properly normalized by the construction, i.e. it must fulfill 

J2 dpik,t,\I,to) = l, (171) 

for an arbitrary starting point (to,/). The probability that the transition to any other 
state occurs before some time t is obtained as 

I Yl dp{K,t,\I,to) = l-e-*^(*'*°), (172) 



to<ti<t 



where we have denoted 



t 

^i{t,to) = j dt'J^Pjiit'). (173) 



to 



The probability that such a transition occurs after the time ti = t is just equal to 



f J2 dp{k,t,\I,to) = e-'''^''''^. (174) 

tit 

All the above is a repetition of the very standard description of the Poissonian "decay 
mechanism" with the "time- dependent" transition (decay) constants Pji{t) > 0; it de- 
scribes precisely what we basically do understand as a Markovian process^. 

Let us now imagine a very large ensemble of identical objects, each of them at a given 
time t in one well defined state k, and evolving statistically independently according to the 
transition probability distribution defined above. Let us introduce the population A^/(t) 
of the objects which are at a given time t in the state /. Given the above probabilistic 
transition rule, we may easily calculate the change of the population ANj{t) from the time 
t to the time t + At. The original population A''/(t) is diminished to A''/(t)e~*^^*"'"^*'*''. At 
the same time interval At the population of the state I is also increased by the influx 
from all other states L ^ I by PuNLAt. Altogether we get 

Nj{t + At) = Ar,(t)e-*^(*+^*'*) + J2PiLit)NLit)At, (175) 



^In the literature one may find many different definitions of tlie Markovian process, see e.g. [23-25] . 
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and 



ANj{t) = Njit + At) - Niit) = -Njit) Pjiit)M + PiL{t)NL{t)At , (176) 
or equivalently 

9tiV,(t) = -(^Pj,(t))iV,(t) + ^P,^(t)iV,.(t). (177) 

We therefore conclude that the differential evolution equation ()169|) can only be compat- 
ible with the probabilistic Markovian process if the following property of the transition 
matrix is true 

Pnit)^-J2Pnit). (178) 

This is what we shall always assume to be true in the standard Markovian process. 

The QCD singlet evolution kernels do not fulfill the above condition. Hence, perfect 
Markovianization (with weight =1 events) is not possible in this context and the use of the 
weighted events is mandatory, at least at the internal level of the parton-shower MC. MC 
events can always be turned into weight=l events with the usual rejection methods, but 
at some price. The remedy is to use zPik{z), which fulfill the above condition, instead of 

PjKiz). 

D.2 Iterative solution 

For completness let us write down the iterative solution of the evolution equation 

mi{t) = -RiNj{t) + Y PiK{t)NK{t), Ri = -Pn , (179) 

in the discrete space. The above equation can be easily brought to a homogeneous form 

e-*^(*'*°)9i(e*^(*'*°)iV,(t)) =YPiKit)NK{t), 



K^I 



^ lit, to) = j dti Ri{ti), 

to 

which then can be turned into an integral equation 

t 

Nr{t) = e-*^(*'*°)iV,(to) + J dt, e'^^^*-*^) P;^(ti) iV;,(ti 



:i80) 



*o ^ (181: 



, ^ Pkj, for K^J, 
~ ^ 0, hi K = J, 
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and finally can be solved by means of multiple iteration 



t 

oo n r- c 

iV^(t)=e-*-(*'*°)iV^(to)+5^ n jdUQ{U-tU, 

n=l Ko,...,Kn-i i=l ^ 



(182) 



X e 



i=l 



NkM, 



where K = Kn, for the brevity of the notation. The above series of integrals with 
positively defined integrands (assuming Pjk > 0) can be interpreted in terms of a random 
Markovian process starting at to and continuing until t. 
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